• 专利附录
  • 专利说明
  • 权利要求
  • 类似技术
  • 同族专利
  • 引用文献
  • 法律状态
  • 优先权
    • 专利摘要:
    Method for calculating the effective permeability of a network of discrete fractures in an underground environment The method may include: obtaining (800) attributes of the discrete fracture network, spatial sampling (810) of points on fracture surfaces within a computing cell, allocation (820) ) at each sampled point of at least one attribute of the corresponding fracture, calculating (850) discrete pressure values in the calculation cell for the location of the sampled points by solving partial differential equations of a model of flow, calculating (860) effective permeability values for the calculation cell from the pressure values.
    公开号:FR3046479A1
    申请号:FR1650015
    申请日:2016-01-04
    公开日:2017-07-07
    发明作者:Emmanuel Malvesin；Jean-Pierre Joonnekindt
    申请人:Services Petroliers Schlumberger SA；
    IPC主号:
    专利说明:

    (équation 2) où : div est l'opérateur de divergence, K est le tenseur de perméabilité du réseau de fractures à l'intérieur de la cellule de simulation 300 ; p est le champ de pression à déterminer,
    Vp est le gradient de pression, p est une pression prescrite appliquée sur la frontière de la cellule de simulation 300, laquelle peut prendre différentes valeurs arbitraires, Ω représente un domaine spatial correspondant aux fractures présentes à l'intérieur de la cellule de simulation 300, Γ représente l'intersection de Ω avec les frontières de la cellule de simulation, c.-à-d. Γ représente l'intersection des fractures présentes à l'intérieur de la cellule de simulation 300 avec les frontières de la cellule de simulation ; il est supposé que les conditions de Dirichlet s'appliquent à Γ.
    [0047] L'équation 1 exprime le principe de conservation de la masse et du volume appliqué au fluide incompressible. L'équation 1 forme les équations aux dérivées partielles (PDE) du modèle basé sur l'écoulement.
    [0048] Le tenseur de perméabilité K est représenté par une matrice symétrique 3x3, dont le coefficient matriciel kÿ représente l'effet directionnel de l'ensemble de fractures présent dans la cellule de simulation 300 sur l'écoulement du fluide dans les trois directions de l'espace relativement au système de coordonnées de simulation (x,y,z) dans lequel la cellule de simulation 300 est définie. Le tenseur de perméabilité K peut donc être écrit sous la forme
    (équation 3) [0049] Les équations 1 et 2 ci-dessus sont applicables aux trois directions spatiales d=l, II, III pour déterminer les valeurs de pression dans ces trois directions d=l, II, III, ces trois directions correspondant respectivement aux axes x,y,z de la cellule de simulation 300.
    [0050] Afin de résoudre numériquement l'équation PDE (équation 1), la forme dite "forte" des équations basées sur l'écoulement représentées par l'équation 1 peut être transformée en la forme dite "faible". La “forme faible” repose sur une approximation de la solution de l'équation 1 au moyen d’une fonction test, notée w(x), et est obtenue en multipliant la forme forte par ladite “fonction test” w(x) et en intégrant par parties la “forme forte” sur le domaine spatial Ω, défini par les fractures dans la cellule de simulation. Comme cela est expliqué ci-dessous, la fonction test w(x) peut être choisie arbitrairement puisque les valeurs de pression calculées p, ne dépendent pas des valeurs de la fonction test.
    [0051] La "forme faible" des équations PDE du modèle basé sur l'écoulement peut être écrite comme suit :
    (équation 5) où la constante a est utilisée pour imposer les conditions aux limites de la cellule de simulation 300 ; a peut être choisie arbitrairement de telle sorte que a » 1 et sert à contrôler les conditions aux limites du domaine spatial Ω ; n est le vecteur unitaire normal à 5ΩΓ ; dS est un élément de surface sur une fracture.
    La méthode "sans maillage" [0052] L'équation PDE (équation 5) peut être résolue de différentes façons afin de calculer les valeurs de pression. Dans un ou plusieurs modes de réalisation, la méthode "sans maillage", qui est décrite de manière plus détaillée ci-dessous, repose sur des méthodes de calcul telles que la méthode dite avec maillage de fond, comme la « méthode Galerkine sans éléments » (ci-après la méthode EFG, Element-Free Galerkin) (divulguée par exemple par Belytschko et al., 1994) ou la méthode purement sans maillage comme la « méthode Petrov Galerkine locale sans maillage» (ci-après la méthode MLPG, Meshless Local Petrov Galerkin) (divulguée par exemple par Atluri et Zhu, 1999). Dans un ou plusieurs modes de réalisation, la méthode "sans maillage" permet de transformer la "forme faible" de l'équation 5 en une forme discrète, et permet donc une discrétisation des équations PDE du modèle basé sur l'écoulement. Cette forme discrète des équations basées sur l’écoulement correspond à un système d'équations linéaires à partir duquel des valeurs de pression discrètes du champ de pression peuvent être calculées.
    Echantillonnage spatial [0053] Afin de générer la forme discrète du modèle basé sur l'écoulement, un ensemble de points échantillonnés P, - également appelés ici “points de simulation” ou “nœuds” - est échantillonné spatialement sur les surfaces des fractures du réseau de fractures. Dans au moins un mode de réalisation, chaque fracture est représentée par une forme plane et une épaisseur et la surface de la fracture correspond donc à la surface de la forme plane représentant cette fracture. La surface de la fracture peut par exemple être représentée par un parallélépipède plan (par exemple un rectangle) et les points Pi sont échantillonnés sur cette surface. Dans au moins un mode de réalisation, afin de tenir compte de l'épaisseur de la fracture, la surface représentant la fracture peut être une forme plane dans un plan médian relativement à l'épaisseur de la fracture.
    [0054] Les figures 5A et 5B illustrent plusieurs modes de réalisation d'échantillonnage de points sur une surface de la fracture. Les deux figures représentent une fracture 310 par un rectangle. Les points échantillonnés P, peuvent être soit un ensemble régulier de points (tel qu’illustré par la figure 5A), c.-à-d. généré en utilisant une période d'échantillonnage spatiale fixe dans une ou deux directions dans la surface de la fracture, soit un ensemble de noeuds irréguliers (tel qu'illustré par la figure 5B), par exemple avec des localisations déterminées de façon aléatoire sur la surface de la fracture. Par exemple, comme l'illustre la figure 5A représentant une fracture rectangulaire, les points échantillonnés Pj peuvent être alignés verticalement et horizontalement sur une grille d'échantillonnage rectangulaire.
    [0055] Les lignes d'intersection entre les fractures n'ont pas à être échantillonnées et l'élimination de tout nœud ne figurant pas dans la cellule de simulation 300 peut être effectuée très simplement et rapidement.
    [0056] Dans un ou plusieurs modes de réalisation, pour chaque point échantillonné Pj, les valeurs des attributs de la fracture à partir de laquelle un point échantillonné Pj a été échantillonné seront attribuées à ce point échantillonné Pj et utilisées dans la résolution des équations aux dérivées partielles. En particulier, la valeur Kt ou la perméabilité induite par la fracture et/ou directement la valeur e de l'épaisseur de la fracture sera attribuée au point échantillonné Pj. Mécanisme d'interpolation - fonctions de forme [0057] Dans un ou plusieurs modes de réalisation, la méthode "sans maillage" repose en outre sur un mécanisme d'approximation pour approximer une fonction (par exemple le champ de pression p) au moyen de fonctions de forme et de valeurs discrètes (par exemple des valeurs de pression p,) associées aux points échantillonnés P,. Différents mécanismes d'approximation peuvent être utilisés, par exemple les moindres carrés mobiles (MLS), la partition de l'unité (PU) ou les fonctions de Shepard. Dans la description ci-dessous, les moindres carrés mobiles (MLS) sont utilisés pour illustrer un exemple de mode de réalisation.
    [0058] Dans le mécanisme d'approximation MLS, une fonction de forme est définie relativement à un point de référence P ayant un localisation spatiale donnée, la fonction de forme ayant par exemple une propriété de symétrie relativement à ce point de référence P. Dans un ou plusieurs modes de réalisation, pour chaque point échantillonné P,, une fonction de forme correspondante est définie en attribuant au point de référence P la localisation spatiale du point échantillonné Pjdans l'aire volumique de la cellule de simulation. Une fonction de forme est ainsi associée à chaque nœud (point échantillonné P,). Par exemple, en supposant une fracture représentée par une forme plane, une fonction de forme bidimensionnelle ¢£^(x,y) peut être définie dans le plan de la fracture pour les points échantillonnés Pt à l’intérieur de cette forme plane. L'ensemble de ces fonctions constitue une base de l'espace fonctionnel dans lequel l'approximation discrète de la solution de l'équation PDE est réalisée.
    [0059] Dans les applications pratiques, la fonction de forme est non nulle au voisinage rapproché du nœud associé (point échantillonné Pj) et ce voisinage est appelé domaine d’influence du nœud. En dehors de ce domaine d'influence, la valeur de la fonction de forme est nulle ou décroissante.
    [0060] Différentes fonctions de forme peuvent être utilisées avec un domaine d'influence circulaire, elliptique, rectangulaire ou de n'importe quelle forme fermée régulière pratique. Le domaine d'influence peut être un domaine d'influence bidimensionnel ou un domaine d'influence tridimensionnel. S'il est tridimensionnel, le domaine d'influence peut être une sphère, une ellipse, un cube ou n'importe quelle autre aire volumique simple. Dans au moins un mode de réalisation, quand une fracture est représentée par une forme plane (par exemple par un rectangle ou un parallélépipède plan), un domaine d'influence bidimensionnel est choisi pour les fonctions de forme. Dans au moins un mode de réalisation, quand le réseau DFN est composé de fractures rectangulaires, un échantillonnage régulier de points sur les surfaces des fractures est exécuté et un domaine d'influence rectangulaire est choisi pour les fonctions de forme.
    [0061] Dans l'exemple illustré à la figure 6A, des domaines d'influence Dj sont définis autour de chaque point échantillonné Pi dans le domaine spatial Ω. Dans l'exemple illustré à la figure 6B, des domaines d'influence Ω, sont définis autour dé chaque point échantillonné Pj dans le domaine spatial Ω.
    [0062] Dans un ou plusieurs modes de réalisation, l'épaisseur de la fracture est prise en compte dans la fonction de forme définie pour un point échantillonné P,. En fait, comme les lignes d'intersection entre deux fractures ne sont pas représentées explicitement par les points échantillonnés, et afin d'améliorer la pertinence de la simulation, l'épaisseur de la fracture est prise en compte dans la fonction de forme dans la direction normale à la fracture. Le paramètre d'épaisseur de la fonction de forme peut être proportionnel ou égal à l'épaisseur attribuée à la fracture dans le réseau DFN. Par exemple, la fonction de forme peut être définie mathématiquement comme suit :
    (équation 8) où
    et où OfLS(x,y) est une fonction de forme bidimensionnelle définie dans le plan de la fracture autour du point Pi et e est l'épaisseur de la fracture. La figure 6C illustre un mode de réalisation d'un domaine d'influence tridimensionnel Ωι d'une fonction de forme tridimensionnelle correspondant à l'équation 8.
    [0063] Dans un ou plusieurs modes de réalisation, le mécanisme d'approximation MLS est utilisé et chaque fonction de forme est alors définie suivant une fonction de pondération wMLS :
    [0064] La fonction de pondération bidimensionnelle wMLS peut par exemple être définie sur la base d’une fonction spline ou d’une fonction gaussienne. Dans au moins un mode de réalisation, étant donné un point échantillonné P; de coordonnées (χϊ,γι) dans le plan de la fracture, la fonction de pondération peut être définie comme suit:
    où
    et r est le rayon d'un domaine d'influence circulaire autour du point de référence Pi- [0065] D’autres exemples de fonctions de forme et de fonctions de pondération sont divulgués dans le document intitulé “An Introduction to mesh-free methods and their programming”, de Liu GR et al., Springer, 2005.
    Discrétisation de l'équation PPE
    [0066] Les fonctions de forme sont utilisées pour approximer le champ de pression p(x,y,z) et la fonction test w(x,y,z) à la localisation (x,y,z) dans la cellule de simulation 300. Dans un ou plusieurs modes de réalisation, un premier ensemble de fonctions de forme peut être utilisé pour l'approximation de la pression, et un second ensemble de fonctions de forme, différent du premier ensemble, peut être utilisé pour l'approximation de la fonction test w(x). Comme autre exemple, le même ensemble de fonctions de forme est utilisé pour les deux approximations.
    [0067] Dans au moins un mode de réalisation, le mécanisme d'approximation MLS est utilisé et l'approximation peut être définie comme suit :
    (équation 9a) et
    (équation 9b) où N est le nombre total de points échantillonnés,
    Wj est une valeur définie pour la fonction test au point échantillonné Ph et Pi est la valeur de pression à déterminer pour le point échantillonné Pi.
    [0068] En injectant l'approximation des équations 9a et 9b dans l'équation PDE de l'équation 5, la forme discrète suivante des équations basées sur l'écoulement est dérivée :
    (équation 10a)
    (équation 10b)
    (équation 10c) où les indices i, j varient de 1 à N.
    [0069] Les équations 10a, 10b, 10c ci-dessus peuvent être également écrites sous la forme :
    Pour tout W e RN, WTAP = WTB où
    A est une matrice N*N de coefficients et B un vecteur N*1 de coefficients Bt :
    [0070] Enfin le système d'équations linéaires à résoudre pour calculer les valeurs de pression pt est indépendant du vecteur W et peut être écrit sous la forme : AP = B (équation 11a) où : "
    (équation 11b)
    (équation 11c) [0071] En supposant que le domaine spatial Ω est composé d'une seule fracture, la perméabilité de cette fracture peut être évaluée par la loi cubique suivante :
    (équation 11 d) où e est l'épaisseur de la fracture et u le vecteur unitaire colinéaire avec la direction principale de la fracture. Comme l'illustre la figure 6c, la direction principale d'une fracture rectangulaire 310 est par exemple la direction dans laquelle la fracture est la plus grande, par exemple la direction de la longueur de la fracture quand la longueur est supérieure à la largeur.
    [0072] Dans un ou plusieurs modes de réalisation, le calcul des intégrales des équations 11b et 11c est basé sur la méthode MLPG. La méthode MLPG est une méthode qui n'utilise aucun maillage pour le calcul des intégrales.
    [0073] Dans d'autres modes de réalisation, le calcul des intégrales des équations 11b et 11c est exécuté en utilisant un maillage de fond à l'intérieur de la cellule de simulation 300. Pour un tel maillage de fond les contraintes exposées ci-dessus dans la section de fond d'un maillage représentant les fractures et leurs intersections n'existent pas. En conséquence, le maillage de fond peut être choisi arbitrairement. Des exemples de maillages de fond pouvant être utilisés pour le calcul d'intégrales sont représentés schématiquement par les figures 7A et 7B. Un maillage de fond triangulaire (figure 7A) ou un maillage de fond rectangulaire/quadrangulaire (figure 7B) peut être utilisé.
    [0074] En utilisant un maillage de fond, chaque intégrale des équations 11b et 11c peut être écrite sous forme de somme calculée sur chaque cellule de fond Elt du maillage de fond. Par exemple, comme chaque intégrale est constante dans la direction z (x, y restant fixes), chaque intégrale sur une aire volumique d'une fonction f(x,y,z) peut être calculée comme la somme d'intégrales calculées sur une cellule de fond Elt
    (équation 12)
    Les équations 11a, 11b et 11c et cette équation 12 sont applicables à une seule fracture ainsi qu'à tout le domaine spatial Ω représentant un ensemble de fractures dans la cellule de simulation 300. Ainsi, le calcul d'intégrales peut être exécuté fracture par fracture avant l'addition des résultats obtenus pour chaque fracture individuelle. Les coefficients A,yef Sy sont donc calculés pour une seule fracture, puis les matrices A et B sont calculées en additionnant les résultats obtenus pour chaque fracture individuelle avant de calculer les valeurs de pression p, en inversant Iq matrice A.
    [0075] Le calcul numérique des intégrales peut être basé sur la formule de quadrature courante (intégration de Gauss) et peut utiliser les points dits de Gauss (voir par exemple le document intitulé “The Finite Elément Method”, de Zienkiewicz, Taylor & J. Z. Zhu, Elsevier, 2005). Par exemple, un maillage de fond quadrangulaire est utilisé avec 1, 2*2 ou 3*3 points de Gauss par cellule quadrangulaire du maillage de fond. La forme générale de la formule de quadrature pour intégrer une fonction f(x) à Ng points de Gauss xG. peut alors être calculée pour une cellule de fond par :
    (équation 13) où wf représente un coefficient de pondération spécifié par la méthode d'intégration gaussienne.
    [0076] Un point de Gauss xG. pourrait appartenir au domaine d'influence de plusieurs fonctions de forme ou tomber en dehors de la cellule de simulation considérée 300. Si un point de Gauss tombe en dehors de la cellule de simulation considérée 300, le point de Gauss est ignoré. Comme l'illustre la figure 6D, quand l'épaisseur e attribuée à une fracture est prise en compte dans la fonction de forme d'un point échantillonné (voir par exemple l'équation 8 ci-dessus), un point de Gauss pourrait appartenir au domaine d'influence de deux fonctions de forme of"·5 et OyMiS. En fait, le fait de tenir compte de l'épaisseur de la fracture dans la fonction de forme permet de simuler la propagation de la pression d'une fracture à une autre par les termes croisés φΜί.5.0mls ou |eurs grac|jents vofL5 · VOfLS présents dans l'équation 11b. En conséquence, l'intersection des fractures n'a pas besoin d'être échantillonnée, ce qui donne un processus de simulation simplifié et un temps de calcul réduit.
    Perméabilité effective [0077] Après le calcul du champ de pression p au moyen d'un ensemble de valeurs de pression discrètes p,, la perméabilité effective peut être calculée pour une cellule de simulation 300. Le calcul peut être exécuté conformément à la méthode divulguée dans le document intitulé “Permeability tensor of three-dimensional fractured porous rock and a comparison to trace map prédictions” de P. S. Lang et al. Journal of Geophysical Research, 2014.
    [0078] Au moyen de cette méthode de calcul, la loi de Darcy est utilisée pour relier la vitesse de filtration u et le gradient de pression :
    (équation 14a) ou
    (équation 14b) où : K est le tenseur de perméabilité effective (matrice 3 * 3) et /¾ les coefficients du tenseur de perméabilité effective K (à déterminer), μ est la viscosité du fluide (constante, peut être choisie égale à 1 ),
    Vp le gradient de pression calculé,
    Uj la vitesse de filtration le long de l'axe j= x, y ou z j^· est le gradient de pression le long de xi dans le système de coordonnées local de la fracture, i,j représentant x,y,z dans le système de coordonnées local de la fracture.
    [0079] Dans au moins un mode de réalisation, dans chaque direction d (où d=l, Il ou III) la version sans maillage des valeurs moyennées en volume sur une cellule de simulation 300 pour un flux orienté o < ut > et le gradient de pression < ^ > est calculée par :
    (équation 15a)
    (équation 15b) où
    est la vitesse de filtration moyenne le long de l'axe /= x, y ou z, est le gradient de pression moyen le long de /
    Elt est un élément d'un maillage de fond défini dans la cellule de simulation 300,
    Vceti est le volume de la cellule de simulation 300.
    [0080] Dans au moins un mode de réalisation, une formule de quadrature peut être utilisée pour calculer l'intégrale (voir par exemple le document intitulé “The Finite Elément Method”, de R. L. Taylor et al., Elsevier, 2005).
    [0081] Le tenseur de perméabilité équivalente ky est donc obtenu par discrétisation de l'équation 14b :
    (équation 15a)
    où : (équation 15b)
    Procédé de calcul de la perméabilité effective [0082] La figure 8 montre un organigramme conformément à un ou plusieurs modes de réalisation d'un procédé de calcul. Bien que les divers blocs de l'organigramme soient présentés et décrits séquentiellement, l'homme de métier réalisera qu'au moins certains des blocs peuvent être exécutés dans des ordres différents, peuvent être combinés ou omis, et qu'au moins certains des blocs peuvent être exécutés en parallèle. Dans au moins un mode de réalisation, le procédé de calcul est exécutée par le système informatique 100. Dans au moins un mode de réalisation, le procédé de calcul est exécuté pour chaque cellule de simulation 300 d'une grille de simulation. La cellule de simulation 300 définit une aire volumique dans la formation souterraine. La cellule de simulation 300 est par exemple un parallélépipède rectangle tridimensionnel.
    [0083] Au bloc 800, un modèle de réseau DFN est obtenu. Le modèle de réseau DFN comprend des attributs des fractures. Dans au moins un mode de réalisation, chaque fracture est représentée dans le modèle de réseau DFN par une forme plane (par exemple un parallélépipède ou un rectangle) et une épaisseur.
    [0084] Au bloc 810, pour chaque fracture dans une cellule de simulation 300, des points P, sont échantillonnés sur la surface de la fracture. Dans au moins un mode de réalisation, les points P, sont échantillonnés sur la surface de la forme plane représentant la fracture. Dans au moins un mode de réalisation, les points sont échantillonnés spatialement de manière aléatoire sur une surface d'au moins une fracture. Dans au moins un mode de réalisation, les points sont échantillonnés à intervalles réguliers dans au moins une direction sur la surface d'au moins une fracture.
    [0085] Au bloc 820, au moins un attribut de fracture est attribué à chaque point échantillonné P,. Dans au moins un mode de réalisation, une épaisseur de la fracture correspondante est attribuée à chaque point ou une perméabilité induite par la fracture correspondante.
    [0086] Au bloc 830, une fonction test est obtenue. La fonction test peut être sélectionnée par un opérateur du système informatique 100 ou peut être sélectionnée par un programme d’ordinateur exécuté sur le système informatique 100. La fonction test peut être sélectionnée parmi un ensemble de fonctions test prédéfinies et/ou des paramètres de la fonction test peuvent être définis par un opérateur ou par un programme d’ordinateur.
    [0087] Au bloc 840, une fonction de forme est obtenue. La fonction de forme peut être sélectionnée par un opérateur du système informatique 100 ou peut être sélectionnée par un programme d’ordinateur exécuté sur le système informatique 100. La fonction test peut être sélectionnée parmi un ensemble de fonctions test prédéfinies et/ou des paramètres de la fonction test peuvent être définis par un opérateur ou par un programme d’ordinateur. La fonction de forme est définie relativement à un point de référence P ayant une localisation spatiale donnée. Pour chaque point échantillonné P,, une fonction de forme correspondante est définie en attribuant au point de référence P la localisation spatiale du point échantillonné Pi dans l'aire volumique de la cellule de simulation 300.
    [0088] Dans au moins un mode de réalisation, une fonction de forme tridimensionnelle est sélectionnée. Dans au moins un mode de réalisation, la fonction de forme tridimensionnelle est une combinaison d'une fonction de forme bidimensionnelle définie dans un plan comportant la forme plane représentant la fracture considérée ou dans un plan parallèle à la forme plane représentant la fracture et d'une fonction de forme unidimensionnelle représentant une épaisseur de la fracture dans une direction normal à la surface de la fracture au point échantillonné considéré.
    [0089] Au bloc 850, des valeurs de pression discrètes p, aux localisations correspondant aux points échantillonnés P, sont calculées en résolvant pour les points échantillonnés Ρ, des équations aux dérivées partielles d'un modèle d'écoulement en utilisant l'au moins un attribut attribué. Dans au moins un mode de réalisation, les valeurs de pression discrètes sont calculées comme solution d'un ensemble d'équations linéaires discrètes obtenu à partir des équations aux dérivées partielles au moyen d'une approximation discrète du champ de pression aux points échantillonnés. Dans au moins un mode de réalisation, les systèmes d’équations discrètes représentent une forme discrète d'un modèle d'écoulement pour les points échantillonnés P,. Dans au moins un mode de réalisation, les valeurs de pression pisont calculées en fonction du système d'équations linéaires des équations 11a, 11b et 11c.
    [0090] Dans un ou plusieurs modes de réalisation, l'approximation discrète repose sur une approximation par les moindres carrés mobiles. Dans au moins un mode de réalisation, l'approximation discrète repose sur un ensemble de fonctions de forme tridimensionnelles définies pour les points échantillonnés, chaque dite fonction de forme tridimensionnelle définissant un domaine d'influence d'un point échantillonné au bloc 840.
    [0091] Au bloc 860, la perméabilité effective est calculée pour la cellule de simulation 300 considérée en fonction des valeurs de pression obtenues au bloc 850. Dans au moins un mode de réalisation, la perméabilité effective est calculée suivant les équations 15a et 15b.
    [0092] Dans un ou plusieurs modes de réalisation, une simulation de réservoir de fluide est exécutée sur la base des valeurs de perméabilité effective calculées au 860.
    [0093] Bien que la description précédente ait été décrite ici relativement à des moyens, matériaux et modes de réalisation particuliers, elle ne se limite pas à ceux-ci. Comme autre exemple, les modes de réalisation peuvent être utilisés conjointement avec un système portatif (c.-à-d., téléphone, ordinateur monté au poignet ou sur l'avant-bras, tablette, ou autre dispositif portatif), système transportable (c.-à-d., ordinateur personnel ou système informatique transportable), un système informatique fixe (c.-à-d., un ordinateur de bureau, un serveur, une grappe, ou un système informatique de haute performance), ou dans un réseau (c.-à-d., un système d'informatique en nuage). Ainsi, les modes de réalisation s'étendent à tous les structures, procédés, utilisations, produits informatiques, et compositions fonctionnellement équivalents entrant dans l'étendue des revendications annexées.
    [0094] Bien que quelques exemples de modes de réalisation aient été décrits de manière détaillée ci-dessus, l'homme de métier réalisera facilement que de nombreuses modifications peuvent leur être apportées. En conséquence, toute modification entrera dans l'étendue de la présente description telle que définie dans les revendications suivantes.
    EXTRAPOLATION OF THE EFFECTIVE PERMEABILITY OF A DISCRETE FRACTURE NETWORK
    PRIOR ART
    [0001] Natural fractures in an underground formation can provide information on the storage of fluids, the flow of fluids, etc. For example, an underground formation with natural fractures may contain a reservoir of fluid. For example, naturally fractured reservoirs account for more than half of the remaining oil reserves. The modeling of the pressure and / or the effective permeability due to a fracture network can therefore provide valuable information for the simulation of a reservoir, for example for well operators, drilling service providers, etc.
    Any numerical simulation of a reservoir can be based on a precise estimate of the effective permeability (also called equivalent permeability). This precision is used for a successful extraction of oil from this tank. In standard processing flows, the calculation begins with a model of the discrete fracture network (DFN) from which the actual permeability can then be calculated.
    [0003] Different methods can be used to calculate the effective permeability. For example, Elfeel's "Static and Dynamic Assessment of DFN permeability upscaling". A. et al., Society of Petroleum Engineers, 2012, describes several methods, particularly the Oda method and the flow-based method. The Oda method is based on the geometry of discrete fracture planes and on a statistical analysis based on the number and size of fractures in each simulation cell of a simulation grid. The flow-based method relies on finite element simulations of the flow model separately in each of the three spatial directions to simulate a single-phase fluid flow due to a pressure gradient across the fracture network in a cell. given simulation of a simulation grid.
    The Oda method is fast but it generally overestimates the equivalent permeability and can be imprecise for loosely connected fracture networks. The flow-based method is much more precise but it is very expensive in terms of calculations. For example, in the flow-based method, a triangular mesh is used to represent DFN fractures as well as fracture intersections. Not only must all fractures be meshed, but all intersections of fractures must exactly correspond to the sides of the triangles of the triangular mesh. By way of example, FIG. 3A represents two rectangular fractures 310, 311 and their intersection 320 as well as a triangular mesh fitted to these two fractures 310, 311 and corresponding exactly to their intersection 320. The triangulation process must absolutely comply the requirements of the finite elements (or finite volume) and the same constraints therefore apply to the simulation cell, which implies that the mesh of the fracture network must be cross-checked for each simulation cell of the simulation grid. It turns out that the process of triangulation (also called mesh process) is not only the most complex task but also the most greedy step in terms of computation time.
    In addition, the optimal size of the simulation cell in the Oda method or the flow-based method is also very difficult to determine and its efficient calculation would ideally require an iterative process.
    RESUME
    Embodiments of the disclosure can provide systems, methods, and computer program product for calculating the effective permeability of a discrete fracture network. The method may include obtaining attributes of the discrete fracture network; spatial sampling of points on the fracture surface within a volume of a computing cell; assigning to each sampled point at least one attribute of the corresponding fracture; calculating discrete pressure values at the locations corresponding to the sampled points by solving partial differential equations of a flow model using the assigned attributes; a calculation of effective permeability values induced by the fractures present in the volume of the calculating cell from the discrete pressure values.
    Still other aspects are provided with a computer program product comprising computer executable instructions which, when executed by a processor, cause said processor to execute an effective permeability calculation method according to any what embodiment disclosed here.
    In still other aspects, there is described a computer system comprising one or more processors for processing data; a memory operably coupled to the one or more processors which includes program instructions for causing said one or more processors to execute an effective permeability calculation method according to any embodiment disclosed herein.
    This summary is provided to present a selection of concepts that are described in more detail below in the detailed description. This summary is not intended to identify key or essential functions of the claimed subject matter or to be used as an aid in limiting the scope of the claimed subject matter.
    BRIEF DESCRIPTION OF THE DRAWINGS
    We can better understand the functions and advantages of the embodiments described with reference to the following description and the accompanying drawings.
    FIG. 1 illustrates a general processing flow for the calculation of the effective permeability in a geological environment; Figures 2A-2B illustrate an example of an embodiment of a computer system; Figures 3A-3B illustrate certain aspects of calculation methods; Figure 4 illustrates some aspects of a calculation method; FIGS. 5A-5B illustrate certain aspects of a calculation method; Figures 6A-6D illustrate some aspects of a calculation method; Figures 7A-7B illustrate some aspects of a calculation method.
    FIG. 8 illustrates an example of a flowchart of a calculation method.
    DETAILED DESCRIPTION
    The description below includes the best embodiment for implementing the implantations described. It will be understood as not limiting but rather illustrative of the general principles of the embodiments.
    In a fracture network model (DFN), large and / or large fractures are explicitly modeled in the form of discrete blocks, whereas an implicit fracture model (IFM) is used to statistically represent the residual portion of the distribution of fractures (ie smaller fractures) by gate properties. A DFN network model includes fracture attributes representative of various fracture structural properties. For example, a DFN model can represent fracture distribution (eg fracture intensity, type of distribution), fracture geometry (eg shape, width, length and thickness) and / or fracture orientation (eg example, average dip, average azimuth). For example, natural fractures may be characterized at least in part by orientation or dip, size, and thickness. For example, a natural fracture may be characterized in part by a length / width ratio.
    As schematically illustrated in FIG. 1, a DFN network model can be calculated from data acquired in the subterranean formation by data acquisition tools. The data acquisition tools may include physical tools configured to obtain measurements of a subterranean formation and to detect certain characteristics of the geological structures of the subterranean formation. For example, data traces 101 can be obtained with the data acquisition tools. The data plots can be analyzed to generate extrapolated well logs 102 which in turn can be used to generate a three-dimensional representation 103 of extrapolated parameters (eg, fracture intensity and / or other fracture attributes). For example, the intensity, distribution and orientation of the fractures will be analyzed according to well logs 102. The density P32 (fracture area / unit volume) can be calculated and extrapolated along the wells. Thus, different statistical functions can be applied to determine the attribute density, the direction and the dip of the fracture, calibrated with respect to the different densities found in the wells.
    The model of a network of discrete fractures (DFN) can be generated according to previously calculated attributes 103 (density P32, orientation of fractures, etc.). Different geometries concerning fracture length distribution can be tested and calibrated with respect to fracture distribution. For example, several fracture attributes can be calculated to generate the DFN model: shape, width, length, thickness, orientation, average dip, average azimuth, etc.
    Finally, properties 104 of extrapolated fractures can be determined for the DFN network model, for example pressure values and effective permeability values.
    [0024] A natural fracture may be a crack or a fracture surface within the rock unrelated to foliation or cleavage in the metamorphic rock along which there has been no movement. A fracture along which there has been a displacement is a fault. When the walls of a fracture have moved only perpendicular to each other, the fracture is called a joint. Fractures can greatly enhance the permeability of rocks by connecting the pores, and for this reason, in some reservoirs, in order to enhance the flow of hydrocarbons fractures are mechanically induced. Fractures can also be called natural fractures to differentiate them from fractures induced by reservoir stimulation or drilling. In some shale reservoirs, natural fractures enhance production by enhancing effective permeability.
    Natural fractures may exist in groups, for example spaced several hundred meters apart in a general direction. Such natural fractures can locally enhance permeability and be advantageous in techniques seeking to improve extraction.
    For example, a naturally fractured reservoir may include a rock matrix and a set of natural fractures. When multiple natural fractures have spread in a formation, multiple natural fractures can form natural fracture networks, which for example may contribute to storage (eg porosity) and fluid flow (eg permeability, transmissibility, etc.). With respect to the production of fluid from such a reservoir, natural fractures may permit comparatively rapid fluid flow and may be present at various scales.
    By way of example, the modeling of the properties of natural fractures can facilitate decision-making by determining the advantageous or harmful nature of the natural fracture for one or more particular purposes. For example, if a natural fracture stores a certain amount of a desired fluid (ie a considerable amount), modeling can help decide where a production well and an injection well can be positioned to extract at least a portion of the desired fluid from the natural fracture (eg, considering the modeled flow due to applied pressure, breakthrough, extraction of the desired fluid, etc.).
    For practical purposes of the dynamic modeling of a reservoir, a simulation grid is defined and the properties of the fractures of the DFN network are extrapolated for each simulation cell of the simulation grid. For example, pressure values and / or actual permeability values can be calculated for each simulation cell. A digital reservoir simulation can then be performed based on the calculated effective permeability values.
    Instead of the triangular mesh of the flow-based method for the simulation of the fluid flow model, the calculation method described here is based on a discrete representation of the fractures inside the volume of a cell. simulation. This discrete representation can be generated by spatially sampling points P, on the surface of the fractures and by assigning to these points P, sampled the attributes of the fracture on which these points Pj were sampled.
    As illustrated in FIG. 3B, and by comparison with the example described previously with reference to FIG. 3A, the sampling of the points P can be carried out in various ways on the surface of the fractures 310, 311. In particular the sampling points P do not need to be aligned along the fracture intersections 320 and can be arbitrarily located on the fracture surfaces.
    The calculation method described here produces a result that is comparatively as accurate as that of the flow-based method but with a considerably reduced computation time, in a ratio of 10 to 100 depending on the complexity of the fracture network. . Pressure values at the spatial locations of these sampled points P, can be calculated by solving a set of linear equations obtained by transforming the partial differential equations of the flow model into a discretized form using approximation functions defined relative to the P points, sampled. This computational method for solving partial differential equations (PDEs) of a flow model and performing discretization is referred to herein as a "non-mesh" method.
    The description presented below refers to illustrative functions, machines, block diagrams and flow charts of the methods, systems and computer program according to one or more exemplary embodiments.
    Each function, machine, block diagram and illustrative flowchart described may be implemented as hardware, software, firmware, middleware, microcode, or any suitable combination thereof. In the case of a software implementation, the functions, machines, blocks of the block diagrams and / or flowcharts can be implemented by computer program instructions or software code, which can be stored or transmitted. on a computer-readable medium, or loaded on a general purpose computer, a specific purpose computer or other programmable data processing apparatus for producing a machine, such that the computer program instructions or software code which are executed on the computer or other programmable data processing apparatus create the means to implement the functions described here.
    Embodiments of computer readable media include, but are not limited to, both computer storage media and communication media having any media that facilitates the transfer of a computer program. from one place to another. Specifically, software instructions or computer readable program code for performing the embodiments described herein may be stored, either temporarily or permanently, in whole or in part, on a non-transitory, computer readable medium of a device. local or remote storage with one or more storage media.
    As used herein, a computer storage medium can be any physical medium that can be read, written or accessed more generally by a computer. Examples of computer storage medium include, but are not limited to, flash drives or other Flash memory devices (eg memory sticks, memory cards, key drives), CD-ROMs or other optical storage means , DVDs, magnetic disk storage devices or other magnetic storage devices, memory chips, RAMs, ROMs, EEPROM memories, smart cards, or any other suitable medium that can be used to carry or store program code in the form of instructions or data that can be read by a computer processor. In addition, various forms of computer readable media may be used to transmit or route instructions to a computer, including a router, gateway, server, or other transmission device, wired (coaxial cable, fiber, twisted pair) , DSL cable) or wireless (infrared, radio, cellular, microwave). The instructions may include code in any computer programming language including, but not limited to, assembly language, C, C ++, Basic, HTML, PHP, Java, Javascript, etc.
    The computer system 100 can be implemented as a single hardware device, for example desktop personal computer (PC), laptop, PDA, mobile phone or can be implemented on devices separate interconnected devices connected to each other by a communication link, with wired and / or wireless segments.
    In one or more embodiments, the computer system 100 operates under the control of an operating system and executes or relies on various computer software applications, components, programs, objects, modules, data structures, etc. .
    As schematically illustrated in FIG. 2A, the computer system 100 comprises a processing unit 110, a memory 111, one or more computer storage media 112, and other associated equipment such as input interfaces. / output (eg device interfaces such as USB interfaces, etc., network interfaces such as Ethernet interfaces, etc.) and a media reader 113 for reading and writing on one or more computer storage media.
    The memory 111 can be a random access memory (RAM), a cache memory, a non-volatile memory, a backup memory (for example programmable memories or flash), read-only memories or any combination thereof. . The processing unit 110 may be any suitable microprocessor, integrated circuit, or CPU with at least one computer-based processing processor.
    In one or more embodiments, the computer storage medium (s) 112 include computer program instructions which, when executed by the computer system 100, cause the computer system 100 to execute one or more of the described methods. here the computer system 100. The processing unit 110 is a hardware processor that processes the instructions. For example, the processing unit 110 may be an integrated circuit for processing the instructions. For example, the processing unit may be one or more cores or microcines of a processor. The processing unit 110 of the computer system 100 may be configured to solicit said one or more computer storage media 112 to store, read, and / or load computer program instructions or software code which, when executed by the processor, causes the processor to execute blocks of a method described herein for the computer system 100. The processing unit 110 may be configured to use the memory 111 when executing the blocks of a method described herein for the computer system 100, for example to load the instructions of the computer program and to store the data generated during the execution of the instructions of the computer program.
    In one or more embodiments, the computer system 100 receives a number of inputs and outputs to communicate information with the outside. As an interface with a user or operator, the computer system 100 generally includes a user interface 114 that incorporates one or more user input / output devices, such as a touch screen, keyboard, mouse, microphone, touchpad, an electronic stylus, or any other type of device. User input can also be received, eg. eg, on a network interface coupled to a communication network, from one or more external computing devices or systems.
    The computer system 100 of FIG. 2A can be connected to a network or be part of it. For example, as shown in FIG. 2B, the network 120 may include multiple nodes (for example the nodes X 122, Y 124). Each node may correspond to a computer system, such as the computer system shown in Figure 2A, or a group of combined nodes may correspond to the computer system shown in Figure 2A. By way of example, the embodiments may be implemented on a node of a distributed system that is connected to other nodes. As another example, the embodiments may be implemented on a distributed computing system having multiple nodes, each portion of one or more embodiments may be located on a different node in the distributed computer system. In addition, one or more elements of the aforementioned computer system 100 may be remotely located and connected to the other network elements.
    Although this is not shown in FIG. 2B, the node may correspond to a blade in a server chassis connected to the other nodes by a backplane. As another example, the node can be a server in a data center. As another example, the node may correspond to a computer processor or microcode of a memory-based computer processor and / or shared resources.
    The nodes (e.g., X-nodes 122, Y-node 124) in the network 120 may be configured to provide services to a client device 126. For example, the nodes may be part of a cloud computing system. The nodes may include functionality to receive requests from the client device 126 and transmit responses to the client device 126. The client device 126 may be a computer system, such as the computer system shown in Figure 2A. In addition, the client device 126 may include and / or execute all or part of one or more embodiments.
    The flow-based model In each simulation cell of a given simulation grid, the effective permeability - also called equivalent permeability - induced by the fractures present in a simulation cell 300 can be calculated. In one or more embodiments, the effective permeability is calculated as a function of pressure values of the pressure field in the simulation cell 300.
    In one or more embodiments, the simulation grid is composed of a set of parallelepiped simulation cells. Each simulation cell 300 is for example a three-dimensional parallelepiped, for example a rectangular parallelepiped. As schematically illustrated in FIG. 4, the attributes of a DFN fracture network are obtained for the fractures present inside a simulation cell 300.
    In one or more embodiments, the calculation of the pressure values and the effective permeability values in a simulation cell of a simulation grid is based on a flow-based model and the equations based on the corresponding flow. These flow-based equations can for example be written by applying the Laplace equations and Darcy's law. To simplify the simulation, the flow-based model is defined for steady-state single-phase flow of an incompressible fluid of constant viscosity through a non-porous medium. In at least one embodiment, the partial derivative equations of the flow-based model can be written as (equation 1)
    (equation 2) where: div is the divergence operator, K is the permeability tensor of the fracture network inside the simulation cell 300; p is the pressure field to be determined,
    Vp is the pressure gradient, p is a prescribed pressure applied on the boundary of the simulation cell 300, which can take different arbitrary values, Ω represents a spatial domain corresponding to the fractures present inside the simulation cell 300, Γ represents the intersection of Ω with the boundaries of the simulation cell, i.e. Γ represents the intersection of the fractures present inside the simulation cell 300 with the boundaries of the simulation cell; it is assumed that the Dirichlet conditions apply to Γ.
    Equation 1 expresses the principle of conservation of mass and volume applied to the incompressible fluid. Equation 1 forms the partial differential equations (PDEs) of the flow-based model.
    The permeability tensor K is represented by a symmetrical matrix 3 × 3, whose matrix coefficient kÿ represents the directional effect of the set of fractures present in the simulation cell 300 on the flow of the fluid in the three directions of the space relative to the simulation coordinate system (x, y, z) in which the simulation cell 300 is defined. The permeability tensor K can therefore be written in the form
    (Equation 3) [0049] Equations 1 and 2 above are applicable to the three spatial directions d 1, II, III for determining the pressure values in these three directions d 1, II, III, these three directions corresponding respectively to the x, y, z axes of the simulation cell 300.
    In order to solve numerically the PDE equation (equation 1), the so-called "strong" form of the equations based on the flow represented by equation 1 can be transformed into the so-called "weak" form. The "weak form" is based on an approximation of the solution of equation 1 by means of a test function, denoted w (x), and is obtained by multiplying the strong form by said "test function" w (x) and by integrating in parts the "strong form" on the spatial domain Ω, defined by the fractures in the simulation cell. As explained below, the test function w (x) can be chosen arbitrarily since the calculated pressure values p do not depend on the values of the test function.
    The "weak form" of the PDE equations of the flow-based model can be written as follows:
    (Equation 5) where the constant a is used to impose the boundary conditions of the simulation cell 300; a can be chosen arbitrarily such that a "1 and serves to control the boundary conditions of the spatial domain Ω; n is the normal unit vector at 5Ω Γ; dS is a surface element on a fracture.
    The "No Mesh" Method The PDE equation (Equation 5) can be solved in different ways to calculate the pressure values. In one or more embodiments, the "no-mesh" method, which is described in greater detail below, relies on computational methods such as the so-called bottom-mesh method, such as the "Galerkine method without elements". (hereinafter the EFG method, Element-Free Galerkin) (disclosed for example by Belytschko et al., 1994) or the purely non-meshed method as the "local Petrov Galerkine method without mesh" (hereinafter the MLPG method, Meshless Local Petrov Galerkin) (disclosed for example by Atluri and Zhu, 1999). In one or more embodiments, the "non-mesh" method makes it possible to transform the "weak form" of equation 5 into a discrete form, and thus allows a discretization of the PDE equations of the flow-based model. This discrete form of flow-based equations corresponds to a system of linear equations from which discrete pressure values of the pressure field can be calculated.
    Spatial Sampling In order to generate the discrete form of the flow-based model, a set of sampled points P, - also referred to herein as "simulation points" or "nodes" - is spatially sampled on the surfaces of the network fractures. fractures. In at least one embodiment, each fracture is represented by a planar shape and a thickness and the surface of the fracture thus corresponds to the surface of the planar shape representing this fracture. The surface of the fracture may for example be represented by a plane parallelepiped (for example a rectangle) and the points Pi are sampled on this surface. In at least one embodiment, in order to take into account the thickness of the fracture, the surface representing the fracture may be a planar shape in a median plane relative to the thickness of the fracture.
    Figures 5A and 5B illustrate several embodiments of sampling points on a surface of the fracture. The two figures represent a fracture 310 by a rectangle. The sampled points P, can be either a regular set of points (as shown in Figure 5A), i.e. generated using a fixed spatial sampling period in one or two directions in the fracture surface, ie a set of irregular nodes (as shown in Figure 5B), for example with locations randomly determined on the fracture surface. For example, as shown in Figure 5A showing a rectangular fracture, the sampled points Pj can be aligned vertically and horizontally on a rectangular sampling grid.
    The intersection lines between the fractures do not have to be sampled and the elimination of any node not appearing in the simulation cell 300 can be performed very simply and quickly.
    In one or more embodiments, for each sampled point Pj, the values of the attributes of the fracture from which a sampled point Pj has been sampled will be attributed to this sampled point Pj and used in the resolution of the equations. partial derivatives. In particular, the value Kt or fracture-induced permeability and / or directly the value e of the thickness of the fracture will be attributed to the sampled point Pj. Interpolation mechanism - shape functions In one or more embodiments, the "non-mesh" method additionally relies on an approximation mechanism for approximating a function (for example the pressure field p) by means of shape functions and discrete values (eg pressure values p,) associated with the sampled points P ,. Different approximation mechanisms can be used, for example least-squares (MLS), partition of the unit (PU) or Shepard functions. In the description below, moving least squares (MLS) are used to illustrate an exemplary embodiment.
    In the MLS approximation mechanism, a shape function is defined relative to a reference point P having a given spatial location, the shape function having for example a property of symmetry relative to this reference point P. In one or more embodiments, for each sampled point P ,, a corresponding shape function is defined by assigning to the reference point P the spatial location of the sampled point Pj in the area of the simulation cell. A shape function is thus associated with each node (sampled point P,). For example, assuming a fracture represented by a planar shape, a two-dimensional shape function £ (x, y) can be defined in the plane of the fracture for the sampled points Pt within this plane shape. The set of these functions constitutes a basis of the functional space in which the discrete approximation of the solution of the PDE equation is realized.
    In practical applications, the form function is non-zero in the close vicinity of the associated node (sampled point Pj) and this neighborhood is called the node's domain of influence. Outside this influence domain, the value of the shape function is zero or decreasing.
    Different shape functions can be used with a circular, elliptical, rectangular or any other convenient closed form of influence. The domain of influence may be a two-dimensional domain of influence or a three-dimensional domain of influence. If it is three-dimensional, the sphere of influence can be a sphere, an ellipse, a cube, or any other simple volume area. In at least one embodiment, when a fracture is represented by a planar shape (for example by a rectangle or a planar parallelepiped), a two-dimensional influence domain is chosen for the shape functions. In at least one embodiment, when the DFN array is composed of rectangular fractures, regular sampling of points on the fracture surfaces is performed and a rectangular influence domain is selected for the shape functions.
    In the example illustrated in FIG. 6A, domains of influence Dj are defined around each sampled point Pi in the spatial domain Ω. In the example illustrated in FIG. 6B, domains of influence Ω are defined around each sampled point Pj in the spatial domain Ω.
    In one or more embodiments, the thickness of the fracture is taken into account in the shape function defined for a sampled point P 1. In fact, as the intersection lines between two fractures are not explicitly represented by the sampled points, and in order to improve the relevance of the simulation, the thickness of the fracture is taken into account in the shape function in the normal direction to the fracture. The thickness parameter of the shape function can be proportional or equal to the thickness assigned to the fracture in the DFN network. For example, the form function can be mathematically defined as follows:
    (equation 8) where
    and where OfLS (x, y) is a two-dimensional form function defined in the plane of the fracture around the point Pi and e is the thickness of the fracture. FIG. 6C illustrates an embodiment of a three-dimensional influence domain Ωι of a three-dimensional shape function corresponding to equation 8.
    In one or more embodiments, the MLS approximation mechanism is used and each shape function is then defined according to a wMLS weighting function:
    The bidimensional weighting function wMLS can for example be defined on the basis of a spline function or a Gaussian function. In at least one embodiment, given a sampled point P; of coordinates (χϊ, γι) in the plane of the fracture, the weighting function can be defined as follows:
    or
    and r is the radius of a circular influence domain around the reference point [0065] Further examples of shape functions and weighting functions are disclosed in the document entitled "An Introduction to mesh-free methods". and their programming ", by Liu GR et al., Springer, 2005.
    Discretization of the EPP equation
    The shape functions are used to approximate the pressure field p (x, y, z) and the test function w (x, y, z) to the location (x, y, z) in the simulation cell. 300. In one or more embodiments, a first set of shape functions may be used for the pressure approximation, and a second set of shape functions, different from the first set, may be used for the approximation of the pressure. the test function w (x). As another example, the same set of shape functions is used for both approximations.
    In at least one embodiment, the MLS approximation mechanism is used and the approximation can be defined as follows:
    (equation 9a) and
    (Equation 9b) where N is the total number of sampled points,
    Wj is a value defined for the test function at the sampled point Ph and Pi is the pressure value to be determined for the sampled point Pi.
    By injecting the approximation of equations 9a and 9b into the PDE equation of equation 5, the following discrete form of flow-based equations is derived:
    (equation 10a)
    (equation 10b)
    (equation 10c) where the indices i, j vary from 1 to N.
    The equations 10a, 10b, 10c above can also be written in the form:
    For all W e RN, WTAP = WTB where
    A is a matrix N * N of coefficients and B a vector N * 1 of coefficients Bt:
    Finally, the system of linear equations to solve for calculating the pressure values pt is independent of the vector W and can be written in the form: AP = B (equation 11a) where: "
    (equation 11b)
    (equation 11c) [0071] Assuming that the spatial domain Ω is composed of a single fracture, the permeability of this fracture can be evaluated by the following cubic law:
    (Equation 11d) where e is the thickness of the fracture and u is the collinear unit vector with the main direction of the fracture. As illustrated in FIG. 6c, the main direction of a rectangular fracture 310 is, for example, the direction in which the fracture is greatest, for example the direction of the length of the fracture when the length is greater than the width.
    In one or more embodiments, the calculation of the integrals of equations 11b and 11c is based on the MLPG method. The MLPG method is a method that does not use any mesh for calculating integrals.
    In other embodiments, the computation of the integrals of equations 11b and 11c is performed by using a background mesh inside the simulation cell 300. For such a background mesh, the constraints set out below above in the bottom section of a mesh representing fractures and their intersections do not exist. As a result, the background mesh can be chosen arbitrarily. Examples of background meshes that can be used for calculating integrals are shown schematically in FIGS. 7A and 7B. A triangular background grid (FIG. 7A) or a rectangular / quadrangular bottom grid (FIG. 7B) can be used.
    Using a background mesh, each integral of equations 11b and 11c can be written as a sum calculated on each bottom cell Elt of the background mesh. For example, since each integral is constant in the z direction (x, y remaining fixed), each integral on a volume area of a function f (x, y, z) can be computed as the sum of integrals calculated on a Elt bottom cell
    (equation 12)
    Equations 11a, 11b and 11c and this equation 12 are applicable to a single fracture as well as to the entire spatial domain Ω representing a set of fractures in the simulation cell 300. Thus, the calculation of integrals can be performed fracture by fracture before the addition of the results obtained for each individual fracture. The coefficients A, yef Sy are thus calculated for a single fracture, then the matrices A and B are calculated by adding the results obtained for each individual fracture before calculating the pressure values p, by inverting Iq matrix A.
    The numerical calculation of the integrals can be based on the current quadrature formula (Gauss integration) and can use the so-called Gauss points (see for example the document entitled "The Finite Element Method", by Zienkiewicz, Taylor &amp; JZ Zhu, Elsevier, 2005). For example, a quadrangular background grid is used with 1, 2 * 2 or 3 * 3 Gauss points per quadrangular cell of the background grid. The general form of the quadrature formula to integrate a function f (x) to Ng Gauss points xG. can then be calculated for a bottom cell by:
    (equation 13) where wf represents a weighting coefficient specified by the Gaussian integration method.
    A point of Gauss xG. could belong to the domain of influence of several form functions or fall outside the considered simulation cell 300. If a Gauss point falls outside the considered simulation cell 300, the Gauss point is ignored. As illustrated in Figure 6D, when the thickness e attributed to a fracture is taken into account in the shape function of a sampled point (see for example equation 8 above), a Gauss point could belong In fact, the fact of taking into account the thickness of the fracture in the shape function makes it possible to simulate the propagation of the pressure of a fracture. another by the crossed terms φΜί.5.0mls or rms vofL5 · VOfLS in equation 11b, so the intersection of the fractures does not need to be sampled, Simplified simulation and reduced calculation time.
    Effective permeability [0077] After calculating the pressure field p by means of a set of discrete pressure values p ,, the effective permeability can be calculated for a simulation cell 300. The calculation can be performed in accordance with the disclosed method in the document titled "Permeability Tensor of Three-Dimensional Fractured Porous Rocks and a Comparison to Map Predictions" by PS Lang et al. Journal of Geophysical Research, 2014.
    By this calculation method, Darcy's law is used to connect the filtration rate u and the pressure gradient:
    (equation 14a) or
    (Equation 14b) where: K is the effective permeability tensor (matrix 3 * 3) and / ¾ the coefficients of the effective permeability tensor K (to be determined), μ is the viscosity of the fluid (constant, can be chosen equal to 1 )
    Vp the calculated pressure gradient,
    Uj the filtration rate along the axis j = x, y or zj ^ · is the pressure gradient along xi in the local coordinate system of the fracture, i, j representing x, y, z in the local coordinate system of the fracture.
    In at least one embodiment, in each direction d (where d = 1, 11 or 11), the non-mesh version of the values averaged in volume on a simulation cell 300 for an o-oriented flow. <ut> and the pressure gradient <^> is calculated by:
    (equation 15a)
    (equation 15b) where
    is the average filtration rate along the axis / = x, y or z, is the mean pressure gradient along /
    Elt is an element of a background mesh defined in the simulation cell 300,
    Vceti is the volume of the simulation cell 300.
    In at least one embodiment, a quadrature formula may be used to calculate the integral (see for example the document entitled "The Finite Element Method", by RL Taylor et al., Elsevier, 2005).
    The equivalent permeability tensor ky is thus obtained by discretization of equation 14b:
    (equation 15a)
    where: (equation 15b)
    Method for calculating effective permeability [0082] Fig. 8 shows a flowchart in accordance with one or more embodiments of a calculation method. Although the various blocks of the flowchart are presented and described sequentially, one skilled in the art will realize that at least some of the blocks may be executed in different orders, may be combined or omitted, and that at least some of the blocks can be run in parallel. In at least one embodiment, the calculation method is executed by the computer system 100. In at least one embodiment, the calculation method is executed for each simulation cell 300 of a simulation grid. The simulation cell 300 defines a volume area in the subterranean formation. The simulation cell 300 is for example a rectangular parallelepiped three-dimensional.
    At block 800, a DFN network model is obtained. The DFN network model includes fracture attributes. In at least one embodiment, each fracture is represented in the DFN pattern by a planar shape (for example a parallelepiped or a rectangle) and a thickness.
    In block 810, for each fracture in a simulation cell 300, points P are sampled on the surface of the fracture. In at least one embodiment, the points P are sampled on the surface of the planar shape representing the fracture. In at least one embodiment, the dots are spatially sampled randomly over an area of at least one fracture. In at least one embodiment, the dots are sampled at regular intervals in at least one direction on the surface of at least one fracture.
    In block 820, at least one fracture attribute is assigned to each sampled point P 1. In at least one embodiment, a thickness of the corresponding fracture is assigned to each point or a permeability induced by the corresponding fracture.
    In block 830, a test function is obtained. The test function may be selected by an operator of the computer system 100 or may be selected by a computer program executed on the computer system 100. The test function may be selected from a set of predefined test functions and / or parameters of the computer. Test function can be defined by an operator or by a computer program.
    In block 840, a shape function is obtained. The shape function may be selected by an operator of the computer system 100 or may be selected by a computer program executed on the computer system 100. The test function may be selected from a set of predefined test functions and / or the test function can be defined by an operator or by a computer program. The shape function is defined relative to a reference point P having a given spatial location. For each sampled point P ,, a corresponding shape function is defined by assigning to the reference point P the spatial location of the sampled point Pi in the volume area of the simulation cell 300.
    In at least one embodiment, a three-dimensional shape function is selected. In at least one embodiment, the three-dimensional shape function is a combination of a two-dimensional shape function defined in a plane comprising the plane shape representing the considered fracture or in a plane parallel to the plane shape representing the fracture and of a one-dimensional shape function representing a thickness of the fracture in a direction normal to the fracture surface at the sampled point under consideration.
    In block 850, discrete pressure values p, at the locations corresponding to the sampled points P, are calculated by solving for the sampled points Ρ, partial differential equations of a flow model using the at least an attribute assigned. In at least one embodiment, the discrete pressure values are calculated as a solution of a set of discrete linear equations obtained from the partial differential equations by means of a discrete approximation of the pressure field at the sampled points. In at least one embodiment, the discrete equation systems represent a discrete form of a flow pattern for the sampled points P i. In at least one embodiment, the pressure values are calculated according to the system of linear equations of equations 11a, 11b and 11c.
    In one or more embodiments, the discrete approximation is based on a moving least squares approximation. In at least one embodiment, the discrete approximation is based on a set of three-dimensional shape functions defined for the sampled points, each said three-dimensional shape function defining a domain of influence of a sampled point at block 840.
    In block 860, the effective permeability is calculated for the simulation cell 300 considered as a function of the pressure values obtained at block 850. In at least one embodiment, the effective permeability is calculated according to equations 15a and 15b.
    In one or more embodiments, a fluid reservoir simulation is performed based on the effective permeability values calculated at 860.
    Although the foregoing description has been described herein with respect to particular means, materials and embodiments, it is not limited to them. As another example, the embodiments may be used in conjunction with a portable system (ie, telephone, wrist or forearm computer, tablet, or other portable device), a transportable system ( ie, personal computer or transportable computer system), a fixed computer system (ie, a desktop computer, a server, a cluster, or a high performance computer system), or in a network (ie, a cloud computing system). Thus, the embodiments extend to all structures, methods, uses, computer products, and functionally equivalent compositions within the scope of the appended claims.
    Although some exemplary embodiments have been described in detail above, those skilled in the art will readily realize that many modifications can be made to them. Accordingly, any modification will fall within the scope of this specification as defined in the following claims.
    权利要求:
    Claims (12)
    [1" id="c-fr-0001]
    A method of calculating the effective permeability of a discrete fracture network in an underground environment comprising: obtaining (800) attributes of the discrete fracture network; spatial sampling (810) of points on the fracture surface within a volume of a computing cell; assigning (820) at each sampled point at least one attribute of the corresponding fracture; calculating (850) discrete pressure values of the locations corresponding to the sampled points by solving partial differential equations of a flow model using the assigned attributes; a calculation (860) of effective permeability values induced by the fractures present in the volume of the calculating cell from the discrete pressure values.
    [2" id="c-fr-0002]
    The method of claim 1 wherein the points are spatially sampled randomly over an area of at least one fracture.
    [3" id="c-fr-0003]
    The method of claim 1 wherein the points are sampled at regular intervals in at least one direction on the surface of at least one fracture.
    [4" id="c-fr-0004]
    The method of claim 1 wherein said at least one attribute of the corresponding fracture comprises a thickness of a corresponding fracture.
    [5" id="c-fr-0005]
    The method of claim 1 or 4 wherein said at least one attribute of the corresponding fracture comprises a permeability induced by the corresponding fracture.
    [6" id="c-fr-0006]
    The method of claim 1 wherein said discrete pressure values are calculated as a solution of a set of discrete linear equations obtained from partial differential equations by discretely approximating a pressure field at the sampled points.
    [7" id="c-fr-0007]
    The method of claim 6 wherein said discrete approximation is based on a moving least squares approximation.
    [8" id="c-fr-0008]
    The method of claim 6 wherein said discrete approximation is based on a set of three-dimensional shape functions defined for the sampled points, a said three-dimensional shape function defining a domain of influence of a sampled point.
    [9" id="c-fr-0009]
    9. The method of claim 8 wherein for each point sampled on a given fracture, said three-dimensional shape function is a combination of a two-dimensional shape function defined in a plane comprising the surface of said fracture and a function of unidimensional shape representing a thickness of said fracture considered in a direction normal to the surface of said fracture considered at the sampled point considered.
    [10" id="c-fr-0010]
    The method of claim 1 further comprising performing a fluid reservoir simulation based on the actual permeability values.
    [11" id="c-fr-0011]
    A computer system comprising: one or more processors for processing data; a memory operably coupled to the processor or processors, including program instructions for causing said one or more processors to perform an effective permeability calculation method according to any one of claims 1 to 10.
    [12" id="c-fr-0012]
    A computer program product comprising computer executable instructions which, when executed by a processor, cause said processor to execute an effective permeability calculation method according to any one of claims 1 to 10.
    类似技术:
    公开号 | 公开日 | 专利标题
    EP2581767B1|2014-04-30|Method for constructing a mesh of a fractured reservoir with a limited number of nodes in the matrix environment
    EP2530493B1|2015-03-18|Method for building a fracture network mesh from a Voronoi diagram
    Siade et al.2010|Snapshot selection for groundwater model reduction using proper orthogonal decomposition
    Hyman et al.2019|Linking structural and transport properties in three‐dimensional fracture networks
    FR3046479A1|2017-07-07|
    Zhang et al.2018|An iterative local updating ensemble smoother for estimation and uncertainty assessment of hydrologic model parameters with multimodal distributions
    FR3041803A1|2017-03-31|
    AU2017342840A1|2019-01-24|Connectivity based approach for field development optimization
    Roubinet et al.2013|Hybrid modeling of heterogeneous geochemical reactions in fractured porous media
    FR3043227A1|2017-05-05|
    FR3034894A1|2016-10-14|
    AU2013393305B2|2017-06-29|2.75D meshing algorithm
    FR3029664A1|2016-06-10|DEFINITION OF NON-LINEAR PETROFACIES FOR A RESERVOIR SIMULATION MODEL
    AU2013393306B2|2017-08-10|2.5D stadia meshing
    AU2013399054B2|2017-08-31|A geostatistical procedure for simulation of the 3D geometry of a natural fracture network conditioned by well bore observations
    CA2821099A1|2014-01-13|Production process for a geological reservoir based on a stock reservoir by calculating an analytical conditional distribution law of uncertain parameters of the model
    US9715762B2|2017-07-25|3D stadia algorithm for discrete network meshing
    EP2813668A1|2014-12-17|Method for optimising the exploitation of a fluid reservoir by taking into consideration a geological and transitional exchange term between matrix blocks and cracks
    FR3036518A1|2016-11-25|INVERSION FOR CONSTRAINTS
    Miranda et al.2014|Natural fracture characterization in Aptian carbonates, Araripe Basin, NE Brazil
    EP3017430B1|2021-04-07|Lofting algorithm for discrete network meshing
    Correia et al.2018|An integrated workflow to combine static and dynamic uncertainties in reservoir simulation models
    US20160245950A1|2016-08-25|Using representative elemental volume to determine subset volume in an area of interest earth model
    Jones et al.2015|Improved Multi-Scale Fracture Models Based on Quantitative Analysis of Outcrop Analogues
    Chugunova et al.2017|Reservoir Simulation Strategies Using Multipoint Statistics for Modeling of Fractured Reservoir Properties
    同族专利:
    公开号 | 公开日
    US20190011600A1|2019-01-10|
    EP3400546A1|2018-11-14|
    FR3046479B1|2018-07-20|
    US11199647B2|2021-12-14|
    EP3400546A4|2019-09-11|
    WO2017120088A1|2017-07-13|
    引用文献:
    公开号 | 申请日 | 公开日 | 申请人 | 专利标题
    FR2981475B1|2011-10-12|2013-11-01|IFP Energies Nouvelles|METHOD FOR CONSTRUCTING A MESH OF A FRACTURE RESERVOIR WITH A LIMITED NUMBER OF NODES IN THE MATRIX|
    US9618652B2|2011-11-04|2017-04-11|Schlumberger Technology Corporation|Method of calibrating fracture geometry to microseismic events|
    WO2014116896A1|2013-01-25|2014-07-31|Services Petroliers Schlumberger|Pressure transient testing with sensitivity analysis|
    US10241232B2|2014-02-03|2019-03-26|Halliburton Energy Services, Inc.|Geomechanical and geophysical computational model for oil and gas stimulation and production|
    US10788604B2|2014-06-25|2020-09-29|Schlumberger Technology Corporation|Fracturing and reactivated fracture volumes|
    KR101506321B1|2014-10-24|2015-03-26|한국지질자원연구원|Method for analyzing multi phase and heat flow of fluids in reservoir and Recording media therefor|
    WO2017082862A1|2015-11-09|2017-05-18|Halliburton Energy Services, Inc.|Fracture network fluid flow simulation with junction area modeling|
    FR3045868B1|2015-12-17|2022-02-11|Ifp Energies Now|METHOD FOR CHARACTERIZING AND EXPLOITING AN UNDERGROUND FORMATION COMPRISING A NETWORK OF FRACTURES|US11125912B2|2013-11-25|2021-09-21|Schlumberger Technology Corporation|Geologic feature splitting|
    CN110688767B|2019-10-09|2021-04-06|浙江大学|Method for evaluating comprehensive difference degree of rock mass fracture network model precision|
    CN111199109B|2020-01-10|2021-09-03|浙江大学|Coupling method for dividing rock mass homogeneous region based on box-counting dimension and cluster analysis|
    法律状态:
    2017-01-27| PLFP| Fee payment|Year of fee payment: 2 |
    2017-07-07| PLSC| Publication of the preliminary search report|Effective date: 20170707 |
    2018-01-25| PLFP| Fee payment|Year of fee payment: 3 |
    2019-11-28| PLFP| Fee payment|Year of fee payment: 5 |
    2020-11-25| PLFP| Fee payment|Year of fee payment: 6 |
    2021-11-25| PLFP| Fee payment|Year of fee payment: 7 |
    优先权:
    申请号 | 申请日 | 专利标题
    FR1650015|2016-01-04|
    FR1650015A|FR3046479B1|2016-01-04|2016-01-04|EXTRAPOLATION OF THE EFFECTIVE PERMEABILITY OF A DISCRETE FRACTURE NETWORK|FR1650015A| FR3046479B1|2016-01-04|2016-01-04|EXTRAPOLATION OF THE EFFECTIVE PERMEABILITY OF A DISCRETE FRACTURE NETWORK|
    PCT/US2016/069067| WO2017120088A1|2016-01-04|2016-12-29|Effective permeability upscaling for a discrete fracture network|
    EP16884217.7A| EP3400546A4|2016-01-04|2016-12-29|Effective permeability upscaling for a discrete fracture network|
    US16/065,812| US11199647B2|2016-01-04|2016-12-29|Effective permeability upscaling for a discrete fracture network|
    [返回顶部]